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Abstract: We perform a detailed dynamical analysis of various cosmological scenarios in 
extended (varying-mass) nonlinear massive gravity. Due to the enhanced freedom in choosing 
the involved free functions, this cosmological paradigm allows for a huge variety of solutions 
that can attract the universe at late times, comparing to scalar-field cosmology or usual 
nonlinear massive gravity. Amongst others, it accepts quintessence, phantom, or cosmological- 
constant-like late-time solutions, which moreover can alleviate the coincidence problem. These 
features seem to be general and non-sensitive to the imposed ansantzes and model parameters, 
and thus extended nonlinear massive gravity can be a good candidate for the description of 
nature. 
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1 Introduction 



The idea of adding mass to the graviton is quite old [1], but the straightforward linear ap- 
proach leads to the van Dam, Veltman, Zakharov (vDVZ) discontinuity [2, 3], that is the 
zero- mass limit of the obtained results does not provide the General Relativity results. This 
is due to the fact that not all the extra degrees of freedom, introduced by the graviton mass, 
decouple at the zero- mass limit, since the longitudinal graviton preserves a finite coupling to 
the trace of the energy-momentum tensor. This discontinuity can be removed if one incorpo- 
rates nonlinear terms [4], however it was soon realized that these necessary nonlinear terms 
introduce the Boulware-Deser (BD) ghost degree-of-freedom [5], making the theory unstable. 

However, recently, a specific nonlinear extension of massive gravity was formulated in 
[6, 7], requiring the Boulware-Deser ghost to be systematically removed (see [8] for a review). 
Such a construction is interesting at the theoretical level, since adding mass to a spin-two 
particle is a well-defined problem by itself, however it has an additional motivation, namely it 
is a new class of (Infra- Red) gravity modification hoping to account for inflation and late-time 
acceleration. The theoretical and phenomenological advantages led to a significant amount 
of relevant research [9-70]. 

Despite the successes of nonlinear massive gravity, it was realized that the usual simple 
homogeneous and isotropic cosmological solutions are unstable at the perturbation level [71], 
which led to less symmetric models [72, 73]. However, in [74] a different approach was followed, 
namely to suitably extend the theory allowing for a varying graviton mass, driven by a scalar 
field. This extended (varying-mass) nonlinear massive gravity proves to exhibit interesting 
cosmological behavior, leading the universe to lie at the quintessence or phantom regime, 
experience the phantom-divide crossing [75], or exhibit bouncing and cyclic behavior [76]. 

Since extended (varying-mass) nonlinear massive gravity exhibits interesting phenomeno- 
logical features when applied to cosmology, in the present work we desire to perform a de- 
tailed dynamical analysis of such a scenario. In this way we can bypass the complexities of 
the equations, which prevent any complete analytical treatment, and investigate in a sys- 
tematic way the huge class of possible late-time cosmological behaviors, calculating various 
observable quantities, such as the dark energy density and equation-of-state parameters and 
the deceleration parameter. 

The plan of the work is the following: In section 2 we briefly review the extended nonlinear 
massive gravity and its cosmological paradigm. In section 3 we perform a dynamical analysis 
of both flat and open geometries, and in section 4 we discuss the cosmological implications 
and the physical behavior of the scenario. Finally, section 5 is devoted to the summary of the 
obtained results. 

2 Cosmology in extended nonlinear massive gravity 

In this section we briefly review cosmology in extended nonlinear massive gravity [74, 75]. 
In this gravitational framework the graviton mass is generalized to be varying, driven by a 
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scalar field. The total action is written as 
"Af|„ „ 



d Xt/— g 



P R + V(iP)(U 2 + a 3 U 3 + aM) - ^d^d^ - Wty) 



+ S m , (2.1) 



where M p is the reduced Planck mass, R is the Ricci scalar, ip is the extra canonical scalar 
field with W{ip) its usual potential and V(tjj) an additional potential coupling to the graviton 
potentials, and a 3 and 04 are dimensionless parameters. The graviton potentials write as 

U 2 = K*JC V V] , U 3 = Kffitfi , U 4 = KffilCyC^ , (2.2) 

with 



K = ^v~ sJgwfABdp^d^ 3 , (2.3) 

where we use the notation 

tcIK] = \{^K-WK) . ( 2 - 4 ) 

and similarly for the other antisymmetric expressions. Furthermore, Jab is a fiducial metric, 
and 4> A (x) are the Stiickelberg scalars introduced to restore general covariance [77]. The above 
extended scenario is still free of the the BD ghost [74]. Finally, in order to obtain a realistic 
cosmology in (2.1) we have allowed for the standard matter action S m , minimally-coupled to 
the dynamical metric, corresponding to energy density p m and pressure p m . 

2.1 Flat universe 

In order to extract the cosmological equations we need to consider specific ansatzes for the 
two metrics. For the physical metric we assume a flat Friedmann-Robertson- Walker (FRW) 
form: 

d 2 s = -N{t) 2 dt 2 + a{t) 2 5 ij dx i dx j , (2.5) 

with a(t) the scale factor and N(t) the lapse function, and for the Stiickelberg fields we 
consider 

0° = b(t), f = a ref x l , (2.6) 

with a re f a (constant) reference scale factor. We mention that contrary to standard massive 
gravity, where such a choice for the dynamical metric cannot be accompanied by a simple 
ansatz for the fiducial one [71], in the present extended scenario the extra freedom does allow 
for a simple Minkowski ansatz for the fiducial metric: 

fAB=VAB- (2.7) 

Variation of the action with respect to iV and a gives rise to the Friedmann equations 

3Mj,H 2 = p DE + Pm , (2.8) 
— 2MpH = pde + PDE + Pm + Pm 5 (2.9) 
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where we have defined the Hubble parameter H = a/a, with a = da/(Ndt), and finally we 
set N = 1. In the above expressions we have defined the effective dark energy density and 
pressure, incorporating the extra gravitational terms, as 

Pde = \i> 2 + W{$) + V{iP) (u - 1) [h(u) + h(u)} (2.10) 

pde = \i? - w{$) - vy>)u(u) - v{ip)bh{u) , (2.ii) 



where 



fi(u) = 3 — 2u + Q3 (3 — u) (1 — u) + CK4 (1 — u) 2 
f 2 (u) = 1 - u + a 3 (1 - uf + ^ (1 - uf 
f 3 (u) = 3 - u + a 3 (1 - u) 

f A (u) = - [6(1 - u) + u 2 + a 3 (1 - it) (4 - 2u) + a 4 (1 - uf , (2.12) 



with 



u = ^l. (2.13) 

a 

These satisfy the usual conservation equation 

Pde + 3H( Pde + Pde ) = 0, (2.14) 

and moreover we can define the dark-energy equation-of-state parameter as 

— "Pde , , ,n 

wde = • (2.15) 

Pde 

Variation of the action (2.1) with respect to the scalar field ift provides its evolution 
equation: 

dW dV 

dtp dip 

Additionally, variation of (2.1) with respect to b provides the constraint equation 

V{^)Hh{u) + V{^)f 2 {u) = . (2.17) 

Finally, one must also consider the matter evolution equation p m + 3H(p m + p m ) = 0. In 
the following we assume matter to have a general equation-of-state parameter w m = 7 — 1 = 
Pm/ Pm, where 7 is the barotropic index, focusing on the usual dust case (7 = 1) only when 
necessary. 

The above cosmological application in a flat universe, although it leads to interesting 
phenomenology, it has significant theoretical disadvantages. These arise mainly from the 
constraint equation (2.17), which using (2.12) in general gives [75]: 

m^lda _ V O a ref 



$ + 3H^ + ^ + ^ {(u - 1) [f 3 (u) + h(u)} + 3£>/ 2 (u)} = . (2.16) 



V(ip(t)) = V e ' »/a(«(»)) 



(a — a Te j)\a^a 2 re ^ — (3«3 + 2af)aa Te f + (3 + 3^3 + a^)a 2 } 

(2.18) 
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As we observe this relation severely restricts the allowed coupling-potential V(ip). Addition- 
ally, as we can see the varying graviton square mass V(ip) diverges and changes sign at least 
for one finite scale factor (namely at a re f), independently of the model parameters, and this 
would make the scenario unstable at the perturbation level. Although one can still choose 
a re j at far past (a re f < 10 -9 ) in order to be smaller than the Big Bang nucleosynthesis scale 
factor and not interfere with the standard thermal history of the universe, or at the far future, 
or even "shield" a re f with a cosmological bounce, case in which the universe is always away 
from it [76], such considerations can only cure the problem phenomenologically, since at the 
theoretical level it remains unsolved. Clearly, the scenario of a flat universe has a serious dis- 
advantage and therefore one should try to construct generalizations in which these problems 
are absent. This will be performed in the next subsection, where the addition of curvature 
makes the graviton mass square always positive. 

2.2 Open universe 

Let us now consider an open 1 FRW form for the physical metric [74]: 



<? s = - N {tfdt 2 + a^tfSad^daP - a(t) 2 k2 ^f dx ^ 2 
w w 3 1 + k 2 (5ijx l x3) 




with N(t) the lapse function and a(t) the scale factor, and K < with k = 
Stiickelberg fields we choose for simplicity the forms 

0° = b(t)Jl + k 2 (5 ij x i xi), <f>* = kbtyx* . (2.20) 

Note that in this case there is no need for the introduction of a reference scale factor ct re j 7 
since it has been absorbed in b(t). Lastly, similarly to the flat case for the fiducial we consider 

fAB=VAB- (2.21) 

Variation of the action (2.1) with respect to N and a gives rise to the following Friedmann 
equations 

3M| (h 2 - ^\ = PDE + p m , (2.22) 

-2M 2 P (h + ^ \ = PDE + pee + p m + p m ■ (2.23) 
In the above expressions we have defined the effective dark energy density and pressure as 
PDE = + W{i>) + V{^) (X - 1) [f 3 (X) + MX)} (2.24) 
PDE = \tf - Wty) - V{i>)U{X) - V{^)bMX) , (2.25) 



1 Similarly to usual massive gravity, closed FRW solutions are not possible since the fiducial Minkowski 
metric cannot be foliated by closed slices [15, 74]. 
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but now the relevant functions become 



h(X) = 3 - 2X + a 3 (3 - X) (1 - X) + a A (1 - X) 2 
f 2 (X) = 1 - X + « 3 (i - xf + ^ (1 - X) 3 
/ 3 (X) = 3-X + a 3 (l-X) 

/ 4 (X) = - [6 - 6X + X 2 + a 3 (1 - X) (4 - 2X) + a 4 (1 - X) 2 



(2.26) 



where 




(2.27) 



These verify the usual conservation equation 



p DE + 3H(p DE + Pde) = 0. 



(2.28) 



Variation of (2.1) with respect to the scalar field tp provides its evolution equation: 



i> + mi, + ^ + ^ {(X - 1) [/ 3 (X) + /x(X)] + 36/ 2 (X)} = . 



(2.29) 



Furthermore, variation with respect to b provides the constraint equation 



Vty) [H--)f 1 (X) + VWf 2 (X) = . 



(2.30) 



Finally, we consider also the matter conservation equation p m + 3H(p m + p m ) = 0. 
3 Dynamical analysis 

In order to investigate the cosmological behavior of the scenario of extended nonlinear massive 
gravity we have to perform its dynamical analysis, and thus we have to transform the involved 
cosmological equations into the autonomous form X' = f(X) [78-81], where X is the column 
vector of suitably introduced auxiliary variables, f(X) the corresponding column vector of the 
autonomous equations, and a prime denotes the derivative with respect to In a. The critical 
points X c are extracted through X' = 0, and in order to examine their stability properties we 
expand around X c as X = X c + U, with U the corresponding perturbations of the variables. 
Thus, at the linear perturbation level and for each critical point we find U' = Q • U, where the 
matrix Q contains the coefficients of the perturbation equations. Therefore, the eigenvalues 
of Q determine the type and stability of the specific critical point. 

The scenario at hand, that was presented in the previous section, consists of the equations 
(2.8), (2.9) or (2.16) and (2.18) for the fiat geometry, and (2.22), (2.23) or (2.29) and (2.30) 
for the open geometry, with a 3 , 04 the model parameters. Although, as we discussed, the 
flat case has theoretical disadvantages, for completeness in the following we analyze it too, 
since it could still be cosmologically valid in suitable frameworks, for example embedded into 
bouncing evolution. 
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As we can see, there are three unknown functions involved, namely the usual scalar 
potential W(ip), the varying graviton mass square V{ip) and the Stiickelberg-field function 
b(t). However, due to the constraint equation, only two out of these three functions are free 
and can be considered as input, while the third one is extracted from the equations of motion. 
As usual, W(tp) is the one function that is always imposed by hand. Throughout the work we 
will consider the usual scalar field potential to have the well-studied exponential form [78-81] 



Thus, in the scenario at hand one could additionally either impose V(ip) at will and leave 
b(t) to be determined by the equations of motion in order to obtain a consistent solution, or 
impose b(t) as an input and leave V(ijj) to be determined by the equations. Definitely, the 
first approach is theoretically more robust, corresponding to the usual Lagrangian description 
where the potentials are imposed as inputs in the theory, and it is the one that is followed 
in all the works on the subject, that is the Stiickelberg fields are always extracted by the 
equations [70-74]. Therefore, in the following subsection we will perform the phase-space 
analysis imposing V(i/j) as an input. However, for completeness, in a separate subsection we 
will also present the (theoretically less interesting) case where b(t) is considered as an input. 

3.1 Imposing V(ip) at will 

For the graviton mass square, and in order to be phenomenologically consistent, without loss 
of generality we assume an exponential form 



In this case the graviton mass is small (at the order of the current Hubble parameter in order 
to drive the current acceleration [8]) at late times, as required by observations, while it could 
play a significant role in the early universe. Additionally, note that in the special case where 
\y = 0, the scenario at hand in the open case corresponds to the usual (constant-mass) 
nonlinear massive gravity. 

3.1.1 Flat universe 

In order to transform the cosmological system (2.8), (2.9) or (2.16) and (2.18) into its au- 
tonomous form, we introduce the dimensionless variables 

aref v W(^) Vty) 

Taking the derivatives of (3.3) and using (2.8), (2.9) and (2.18), we obtain the evolution 
equations for u, Y, and Z, that is the autonomous form of the cosmological system, as 



W(ip) = W e 



(3.1) 



-\vi> 



(3.2) 



u = —u 




(3.4) 
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where primes denote derivative with respect to In a. In the above expressions q = —1 — 
is the deceleration parameter, and the involved H can be expressed in terms of the auxiliary 
variables as 

■ _ H 2 9l (u,Y,Z) 

with 

9l (u,Y,Z) = uh - /2^) + 3A|rZ/f {3/ 2 / 4 - (it- l)(/i + / 3 )[/i - 3( 7 - l)/ 2 ]} 

+3/2/2 + 3Ay Y/f (3 7 A y / 2 - X w h) - 9 7 A 2 y /| + ( 7 - 2)^f 2 , (3.6) 

as it arises from (2.9) and (2.16) through elimination of b (for simplicity we have omitted the 
argument u in fi(u) and f 2 (u)). On the other hand, H elimination between (2.9) and (2.16) 
gives 

g 2 {u,Y,Z) 

with 



^= .g v r'%, . (3-7) 



52 (n, y, Z) = "2n/ 2 ^ + /i f 2u^ - 3 7 / 2 J + 6/i/ 2 + 3y/ 2 ( 7 A - 2\ v X w f 2 ) 

+3Zf 2 {(u - l)(/x + / 3 ) [( 7 - 1)/! - 2X 2 v f 2 ] + hh) + (3.8) 

In summary, (3.4) accounts for an autonomous system defined in the phase space 

{u, Y,Z):0< £^* 2 + (u- l)Z[h{u) + h{u)] + Y < l,u > 0, Y > 0, Z > o) , (3.9) 
6Af/ 2 (u) 2 J 

which is in general non-compact. Furthermore, using (2.8) we can express the dark energy 
density parameter Qde = i n terms of the auxiliary variables as 

while using (3.5) and (3.7) we can express the dark energy equation-of-state parameter and 
the deceleration parameter respectively as 

7 f / 1 (u) fla (u,y,z) f , A \ h(u) 2 y 

WDE = TT-vi ; (3-H) 

(X - l)Z[h(u) + / s (u)] + + Y 

§^ J, - (3-12) 



6A^ - #/ 2 
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Or. sr. 


u c 






Existence 


Stable 101 




"BE 


W DE 


q 


Pi 











for Ay > | 


7<min{l,^},A y >| 
saddle point otherwise 


3 

2^ 


7-1 


22-1 
2 1 


P2 








3-2A^, 


ju > 0,0 < A^ < | 
or /i < 0, A y > | 


7 >1,M>1 
saddle point otherwise 


1 





1 
2 


p 3 





2AJ- 





for A v > | 


^<min{l )7 } ! A y >| 
saddle point otherwise 


1 


Ayy -I 
Ay, L 


3 Aw I 
2A V - - 1 



Table 1. The real and physically meaningful critical points of the autonomous system (3.4), their 
existence and stability conditions, and the corresponding values of the dark-energy density parameter 
Qde, of the dark-energy equation-of-state parameter wde, and of the deceleration parameter q. We 
have introduced the notation fi = (4«3 + + 6). 



The real and physically meaningful critical points (u c , Y C ,Z C ) of the autonomous system 
(3.4) (that is corresponding to < O/xe < 1), are obtained by setting the left-hand-sides of 
these equations to zero, and they are presented in Table 1, along with their existence condi- 
tions. For each critical point we calculate the 3x3 matrix Q of the linearized perturbation 
equations of the system (3.4), and examining the sign of the real part of the eigenvalues of Q 
we determine the type and stability of this point. The details of the analysis and the various 
eigenvalues are presented in Appendix A.l, and in Table 1 we summarize the stability results 
(note that in the case of standard matter (7 = 1) points P\ and P2 belong to a curve of 
critical points). Finally, using (3.10), (3.11) and (3.12), for each critical point we calculate 
the corresponding values of £Ide, wde, and q. 

3.1.2 Open universe 

In order to transform the cosmological system (2.22), (2.23) or (2.29) and (2.30) into its 
autonomous form, we introduce the dimensionless variables 

x- kb Y - ww z- vw 11- * n„- k rtl^ 
T' W z ~ ajj5-' u ~^' Uk " -oh- (3 - 13) 

Differentiating with respect to In a we obtain the autonomous form of the cosmological system: 
X' = -X + Q k b 



Y' = Y 



2(9 + 1) - VE\ w u 



Z' = Z 2(q + 1) - \/6Ay£/] 

U' = 3^X v Zf 2 b + i [V6 [X v (X - l)Z{h + h) + X W Y] + 2(q - 2)u} 

n' k = g n k , (3.14) 
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(3.16) 



with q = — 1 — jp, and where for simplicity we have omitted the argument X in fi(X) and 
f 2 {X). In the above expressions H and b are given by (2.23) and (2.29) as 

H 2 gi (X,Y,Z,U,Q k ,H 2 ) 
H= wu ' ' ' ' fe ' '- r (3.15) 

Xv |_2(o fc - \)n k f 2 $k + 2(n k - i)n fc /ig + 3Zf 2 f 2 [6(n k - i) 2 h 2 + 1] j 

h= g 2 (X,Y,Z,U,n k ,H 2 ) 

Xv {-2(n k - l)n k f 2 ^ + 2(n k - l)f2 fc /ig + 3Z/2/ 2 [6(fi fc - 1)2^2 + !] j ' 

with 

ffl (x, y, z, u, n k ,H 2 ) = J| {-3\ v (n k - \)zhh \{i -i)(x- i)n k + x] 

-Ay(O fc - l)fi fc / 2 {3Z [( 7 - - l)/ 3 + / 4 ] + 3( 7 - 2)U 2 + 3 7 (Y + ^ - l) - 2ft 2 }} 
+|| {3Ay(O fc - 1)Z/ 2 [(7 - - l)n k + X] 

-Av(fi fc - l)Q k fi [-3(7 - 1)(X - l)Zh ~ ZZU ~ 3(7 - 2)U 2 - 3 7 (y + Q 2 k - l) + 2fijJ] } 
+/1V2 {9Ay(0 fe - 1) 2 ZH 2 {3Z [(7 - - l)/ 3 + U] + 3( 7 - 2)U 2 

+3 7 (y + n 2 k - i) - 2n 2 k } - z\ v n k z} 

+ff {27( 7 - l)Xv(X ~ - l) 2 Z 2 / 2 tf 2 

+9(ft fc - 1) 2 Z# 2 -A F (X - l)Z/ 3 + \/6C/ - \ W Y } 

-9A y (X - l)(n fc - l) 2 Z 2 f?H 2 (3.17) 

<? 2 (x, y, z, u, n k ,H 2 ) = -2\ v x(n k - 1)/ 2 ^ + 2A y x(ft fc - i)h^ 

dX aX 
-Av/1/2 {3Z [( 7 - 1)(X - l)/ 3 + / 4 ] + 3( 7 - 2)C/ 2 + 3 7 (Y + ft 2 , - l) - 2(ft fc - l)ft fe } 

+/ 2 {e(ft fc - i) 2 tf 2 [-x v (x - l)Zh + V6U- \ w y\ - 3( 7 - l)\ v {X - l)Zf 2 } 

-6X V (X - l)(ft fe - l) 2 ZffH 2 , (3.18) 

where H 2 is given from (2.30) as 



H 2 



(3.19) 



.V6(l-^)/i(X)C/. 

In summary, (3.14) accounts for an autonomous system defined in the phase space 

{(X,Y,Z,U,n k ) :0< (X-l)Z[h(X) + f 3 (X)] + U 2 + Y + n 2 k < 1,X>0,Y >0,Z>0}, 

(3.20) 

which is in general non-compact. Furthermore, using (2.22) we can express the dark energy 
density parameter in terms of the auxiliary variables as 

V DE = |g§ = (X - l)Z \h{X) + h{X)\ + U 2 + Y, (3.21) 
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Exists for 


Qx 





Ay 


Z c 










Ay— Aiy 


< - ^Zc < 1 


g 2 








Z c 


Vi 

Av 





< 2Xj " ^ < 1 


Q 3 








Z c 


V6 





< - iiZ c < 1 


Qi 











f/c 





< (7 2 < 1 


Q 5 





y c5 










n ^ 3X v +6X v U 2 ~VEU c (\ v X w +3) ^ 1 
U - 3(X V -Xw) ~ 1 


Qt 


x c 


4 





3A[y 




A^>2 


Qt 








4 

3A^/j 


V6 
3A,- 




A?, > 2 


Qs 


1 











1 


always 


Qg 


2az -\-0L4 — yj 4a| — 6a4 

Q4 











1 


a 2 > |a 4 


Qw 


2a3+a4+\/ 4a|— 6a4 











1 


a 2 > |o4 


Qn 













-1 





Table 2. The real and physically meaningful curves of critical points, and individual critical points, of 
the autonomous system (3.14) and their existence conditions, for the case of dust matter (7 = 1). We 
have introduced the notations a = (4^3 + 04 + 6), Y c $ = ^ Xv+Uc { 2X vXw+\ / 6UAx v +x w ) 6] ^ _ 
is the unique real solution of the equation — 2a-$ (X 2 + X — 2) + 04 (X + l)(X — l) 2 + 6 = 0, and 

W \ Ji > - [3a 3 {X-l)-a 4 (X-iy 2 -3]{4a 3 +a 4 +Xla 3 (X-5)+ ai (X-2)-3]+6}- 

while using (3.16) we can express the dark energy equation-of-state parameter as 



w D E 



h{X)b + U(X) 



+ U 2 -Y 

(3.22) 



(X - l)Z[h(X) + h{X)\ + U 2 + Y' 
and finally using (3. 15), (3. 19) the deceleration parameter is expressed as 

= _ 1 9l (X,Y,Z,U,n k ,H 2 ) 

Ay {-2(n k - i)n fc / 2 g. + 2(n k - i)n fc /ig + 3Z/f/ 2 [e(n fc - i) 2 # 2 + 1]} ' 

(3.23) 

Let us extract the critical points of the autonomous system (3.14), setting the left-hand- 
sides of these equations to zero. From the last equation of (3.14) it follows that either q = 
or 0, k = 0, and therefore we can simplify the investigation and examine these two cases 
separately. The details of the analysis, the critical points and critical curves, the various 
eigenvalues and the stability conditions are presented in Appendix A. 2, and in the Table 2 we 
display the real and physically meaningful critical points and their existence conditions for 
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Cr. P. 
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Non-Hyperbolic, 3D stable manifold 
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2 
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Non- Hyperbolic, 3D stable manifold for 

-1 < U c < |fg < A v - < y/%U c , X w < X^{U c ,Xv) or 

1 ^11 ^ n fell T \ / 3 rr \ \ * IT T \ \ 
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Table 3. The stability conditions and the values of the observablcs QdEi wf>e and q, for the 
real and physically meaningful curves of critical points, and individual critical points, of the au- 
tonomous system (3.14), for the case of dust matter (7 = 1). The notations are the same with 
Table 2. Additionally, we have defined w DE i = A v[z e (4 a3 + a4 + 6 A )-7]-1 w z c (4 a 3+a4+6) ' w de 5 = 

3(Av-A t v)[v / 6Av + v / 6A w C/ e 2 -C/ e (2AvA w +3)] 3 A v +6 A v U 2 - V6t7c(AyA t v +3) , (JJ * n 

(3[/-V6Av)[3Av+6AvC/=-V6;7 c (AvA V v+3)] ' U ° E5 ~ 3(Ay-A w ) ' anCl A wK U o^V) - 

V6A 2 , (U^ + l) -3A v (U 2 +3) U a +3VeU^ 
U c (2X 2 v +3uf-2VeX v Uc) ' 



the most interesting case of dust matter (7 = 1), while in Table 3 we present their stability 
conditions and the values of the observables Qde, wde, and q using (3.21), (3.22) and (3.23). 

We mention here, that the variable choice (3.13) allows for an easy, partial, classification 
of expanding and contracting solutions. In particular, solutions with fifc = k/(aH) > 
correspond to H > and thus to expansion, while those with Sl^ < correspond to H < 
and therefore to contraction (k = \/\K\ throughout this work). That is why points with 

> are denoted with the subscript "+", while those with f2fc < are denoted with the 
subscript "-" . However, this is only a partial classification, since it cannot work for solutions 
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Table 4. The interesting individual critical points of the curve of critical points Q5 of Table 3, their 
existence and stability conditions, and the corresponding values of the observables £Ide, wde, and q. 

with = 0, which can be either expanding or contracting. Furthermore, note that although 
our model admits expanding and contracting solutions, from the fifth equation of (3.14) we 
deduce that the sign of Q,k is invariant, and thus transitions from contracting to expanding 
solutions or vice versa do not exist. Nevertheless, since such transitions do exist in the flat 
geometry [76], there could still exist in the non-flat scenario at hand too, but at the edge 
of the phase space, which could be revealed only through application of Poincare central 
projection method [82-84]. This analysis lies beyond the scope of the present work and it is 
left for future investigation. 

Finally, we stress that the curve of critical points Q5 contains many interesting individual 
points, and for that reason we display them separately in Table 4, along with their existence 
and stability conditions and the corresponding values of the observables. Note that these 
points contain the standard quintessence points [79, 85], however the stability conditions are 
slightly different, due to the presence of extra phase-space dimensions, namely curvature and 
graviton mass. 

3.2 Imposing b(t) at will 

In the previous subsection we performed the dynamical analysis following the theoretically 
robust approach in Lagrangian descriptions, that is imposing the potential V(ip) (graviton 
varying square mass) as an input and letting the Stiickelberg field function b(t) to be deter- 
mined by the equations. However, for completeness, and in order to compare with similar 
studies in the literature [86], in this section we follow the theoretically less justified, alterna- 
tive approach, that is to impose b(t) at will and let V(ip) be determined by the equations. 
Similarly to the previous subsection, we will consider the flat and open geometry separately, 
using different b(t) ansantzes in the two case for convenience. 
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3.2.1 Flat universe 

In this case we impose b(t) = Bt with B > 0, since this leads to b = B which simplifies signif- 
icantly the cosmological equations. In the following we focus on the dust matter case (7 = 1), 
however the analysis can be straightforwardly extended to the general 7 case too. In order to 
transform the cosmological system (2.8), (2.9) or (2.16) and (2.18) into its autonomous form, 
we introduce the dimensionless variables 



Taking derivatives with respect to In a, we obtain the autonomous form of the cosmological 
system as 

x' = (q-2)x+^\y 2 

3v [3a 3 + a 4 + u 2 (a 3 + a 4 ) - 2u(2a 3 + a 4 + 1) + 3] {4a 3 + a 4 + u [a 3 (u - 5) + a 4 (u - 2) - 3] + 6} 



2x [-3a 3 (« - 1) + a 4 {u - l) 2 + 3] 
3Bv [3a 3 + a 4 + u 2 (a 3 + a 4 ) - 2u(2a 3 + Q4 + 1) + 3 



2x 



(3.25) 



y' = y(q- \ll Xx + 1 J > ( 3 - 26 ) 

u = -u, (3.27) 

v' = |(u — 1) [— 3a 3 (u — 1) + a 4 (u — l) 2 + 3] 2 | {3v (3a 3 + a 4 + a 4 u 2 — 3a 3 u — 2a 4 u + 3) 

x [3a 3 + a 4 + u 2 (a 3 + a 4 ) - 2u{2a 3 + a 4 + 1) + 3] } + 2(q + l)v. (3.28) 
The autonomous system (3.25)-(3.28) defines a flow in the phase space 

(x,y,u,v) : < - {(u - l)v{u [(u - 5)a 3 + (u- 2)a 4 - 3] +4a 3 + a 4 + 6} + 3 (x 2 + y 2 )} < 1, 

O 

(u — l)v (u 2 a 4 — 3ua 3 — 2ua 4 + 3a 3 + a 4 + 3) I , 

^ = J - <0,u>0,v>0} , 3.29 

where the second inequality follows from the requirement the graviton mass square V(ijj) to 
remain positive. Finally, using (2. 8), (2. 9) and (2.15) we can express the dark energy density 
parameter, the dark energy equation-of-state parameter and the deceleration parameter, in 
terms of the auxiliary variables respectively as 

n DE = -{(u-l)v {u [{u - 5)a 3 + {u- 2)a 4 - 3] + 4a 3 + a 4 + 6} + 3 (x 2 + y 2 )}, (3.30) 

(J 

v [4a 3 + a 4 + n 2 (2a 3 + a 4 + 1) - 2n(3a 3 + a 4 + 3) + 6] + x 2 - y 2 

vjde ~ 



(u — l)v {4a 3 + a 4 + u [a 3 (n — 5) + a 4 (n — 2) — 3] + 6} + x 2 + y 2 
Bv [3a 3 + a 4 + n 2 (a 3 + a 4) — 2n(2a 3 + a 4 + 1) + 3] 
(u — l)v {4a 3 + a 4 + u [a 3 (u — 5) + a 4 (u — 2) — 3] + 6} + x 2 + y 2 



(3.31) 
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Table 5. The real and physically meaningful critical points of the autonomous system (3.25)-(3.28) and 
their existence conditions. We have introduced the notations y,\ = [Aa^ + a 4 — B(3a s + a 4 + 3) + 6] 

nnrl n„ — 403+04+6 

anu ^ 2 - 4a3+04--B(3a3+a 4 +3)+6- 

q = 1 |3w [4a 3 + a 4 + ti 2 (2a 3 + a 4 + 1) - 2u(3a 3 + a 4 + 3) + 6] + 3x 2 - 3y 2 + 1} 
3 

[3a 3 + a 4 + u 2 (a 3 + a 4 ) - 2u(2a 3 + a 4 + 1) + 3] . (3.32) 

The real and physically meaningful critical points (x c , y c ,Uc, v c) of the autonomous system 
(3. 25)- (3. 28), along with their existence conditions, are presented in Table 5. For each critical 
point we calculate the 4x4 matrix Q of the linearized perturbation equations, and we 
determine its type and stability by examining the sign of the real part of the eigenvalues of 
Q. The details of the analysis and the various eigenvalues are presented in Appendix B.l, and 
in Table 6 we display the stability conditions and the corresponding values of the observables 
HflE, wde and q. 

We mention here that the variable choice (3.24) allows for an easy classification of ex- 
panding and contracting solutions. In particular, solutions with y > correspond to H > 
and thus to expansion, while those with y < correspond to H < and therefore to contrac- 
tion. That is why points with y > are denoted with the subscript "+", while those with 
y < are denoted with the subscript However, from (3.26) it is implied that the sign of 
y is invariant, and thus transitions from contracting to expanding solutions or vice versa do 
not exist (there could still exist at the edge of the phase space, which could be revealed only 
through application of Poincare central projection method [82-84], but such an analysis lies 
beyond the scope of the present work and it is left for future investigation). 
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Table 6. The stability conditions and the values of the observables £Ide, wde and q, for the real and 
physically meaningful critical points of the autonomous system (3.25)- (3.28). The notations are the 
same with Tabic 5. 



3.2.2 Open universe 

In this case it proves convenient to impose the ansatz b(t) = boa(t), since this leads to b = boa, 
which simplifies significantly the cosmological equations. In the following, we focus on the 
dust matter (7 = 1), however the analysis can be straightforwardly extended to the general 
7 case too. In order to transform the cosmological system (2.22), (2.23) or (2.29) and (2.30) 
into its autonomous form, we introduce the dimensionless variables 



1/, y/ggj k Vty) k 



VEH ,V V3H a H 2 aH' 

where k = y/\K\. Taking derivatives with respect to In a we obtain the autonomous form of 
the cosmological system as 

, _ 3u 2 vxp5 3u 2 vx5 [2a 3 (P 2 + p - 2) - q 4 (/3 + l)(/3 - l) 2 - 6] n 2 x 3x 3 
n 3 k n 2 [3a 3 (/3-l)-a 4 (/3-l) 2 -3] 2~ + ~Y 

[3 2 . vxS f 6 U 2 [a 3 (/3-4)(/3-l) + a 4 (/3- l) 2 -3/3 + 6] \ 
+ i~2 V Xw+ m\ 3a 3 (/3-l)-a 4 (/3-l) 2 -3 P ] 

+^x {v [/3 2 (2a 3 + a 4 + 1) - 2/3(3a 3 + a 4 + 3) + 4a 3 + a 4 + 6] - 3 (y 2 + l) } , (3.34) 
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y' = \y {v [/? 2 (2a 3 + 04 + 1) - 2/3(3a 3 + « 4 + 3) + 4a 3 + a 4 + 6] - 3y 2 + 3} 

vy(35 nly 3x 2 y [3 x 

7T — + - \ -xy\ w , (3.35) 

2Sl k 2 2 V 2 y ' v ; 

it = —it, (3.36) 

, v 2 f35 ^ 2 18u 2 to 2 <5 

u = — vili. + 3vx + 



n k k ' ft 2 (/3-l)[-3a 3 (/3-l) + a 4 (/3-l) 2 + 3] 

+u [/3 2 (2a 3 + a 4 + 1) - 2/3(3a 3 + a 4 + 3) + 4a 3 + a 4 + 6] - 3y 2 + 3} 
I8u 2 vx 2 5 



(3.37) 



Q k (f3 - 1) [-3a 3 (/3 - 1) + a 4 (/3 - l) 2 + 3] ' 

f 1 3x 2 1 

fi'fc = ^fe | 2 [ /3 ^ 2a3 + «4 + 1) - 2/3(3a 3 + a 4 + 3) + 4a 3 + a 4 + 6] - 3y 2 + l} + — \ 

- V J^-% (3.38) 
2 2 v ' 

where /3 = b^k and 5 = /3 2 (a 3 + q 4 ) — 2/3(2a 3 + a 4 + 1) + 3a 3 + a 4 + 3. 
The autonomous system (3.34)-(3.38) defines a flow in the phase space 

(x, y, u, v, n k ) : < 1 {«(/3 - 1) [a 3 (/3 - 4)(/3 - 1) + a 4 (/3 - l) 2 -3/3 + 6] 

+3{x 2 + y 2 )} +n 2 k < l,u > 0,« > 0} . (3.39) 

Finally, using (2. 22), (2. 23) we can express the dark energy density parameter, the dark energy 
equation-of-state parameter and the deceleration parameter, in terms of the auxiliary variables 
respectively as 

Vde = \ {v(P - 1) [a 3 (0 - 4)(0 - 1) + (u{fi - l) 2 - 3/3 + 6] + 3 {x 2 + y 2 ) } , 

- [/3 2 (2a 3 + q 4 + 1) - 2/3(3a 3 + a 4 + 3) + 4a 3 + a 4 + 6] + 3 (x 2 - y 2 ) 
(/3 - 1> [a 3 (/3 - 4)(/3 - 1) + a 4 (/3 - l) 2 - 3/3 + 6] + 3 (x 2 + y 2 ) 

fiv [/3 2 (a 3 + a 4 ) - 2/3(2a 3 + a 4 + 1) + 3a 3 + a 4 + 3] 

" n k {(/3 - l)v [a 3 (/3 - 4)(/3 - 1) + a 4 (/3 - l) 2 - 3/3 + 6] + 3 (x 2 + y 2 )} ' 
_ /ft; [/3 2 (a 3 + a 4 ) - 2/3(2a 3 + a 4 + 1) + 3a 3 + a 4 + 3] _ n| 

20^. 2 

+- {t> [/3 2 (2q 3 + a 4 + 1) - 2/3(3a 3 + a 4 + 3) + 4a 3 + a 4 + 6] + 3x 2 - 3y 2 + 1} . (3.40) 

The real and physically meaningful critical points (x c ,y c ,u c ,v c ,{l kc ) of the autonomous 
system (3.34)-(3.38), along with their existence conditions, are displayed in Table 7. For each 
critical point we calculate the 5x5 matrix Q of the linearized perturbation equations, and we 
determine its type and stability by examining the sign of the real part of the eigenvalues of Q. 
The details of the analysis and the various eigenvalues are presented in Appendix B.2, and 
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Table 7. The real and physically meaningful critical points of the autonomous system (3.34)-(3.38) 
and their existence conditions. 
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Table 8. The stability conditions and the values of the observables £Ide, wde and q, for the real and 
physically meaningful critical points of the autonomous system (3.34)-(3.38). 

in Table 8 we display the stability conditions and the corresponding values of the observables 
^DE, wde and q. 

Note that the variable choice (3.33) allows for an easy classification of expanding and 
contracting solutions. In particular, solutions with 0/% = k/(aH) > or y > correspond to 
H > and thus to expansion, while those with < or y < correspond to H < and 
therefore to contraction. However, from the equations (3.35) and (3.38) we deduce that the 
sign of y and the sign of £lk are invariant and thus transitions from contracting to expanding 
solutions or vice versa do not exist. Nevertheless, since such transitions do exist in the flat 
geometry [76], there could still exist in the non-flat scenario too, at the edge of the phase 
space, which could be revealed only through application of Poincare central projection method 
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[82-84]. This analysis lies beyond the scope of the present work and it is left for a future 
project. 

4 Cosmological Implications 

In the previous section we performed a complete dynamical analysis of the scenario of extended 
(varying mass) nonlinear massive gravity for both flat and open FRW geometries, we extracted 
the late-time stable solutions and we calculated the corresponding observables. In this section 
we discuss the cosmological implications of the various scenarios case by case. 

4.1 Imposing V(ip) at will 
4.1.1 Flat universe 

First of all we mention that the scenario at hand coincides with standard quintessence if 
the graviton mass square V(i(j) = Voe~ Xv ^ is identically zero. If this is not the case then 
standard quintessence can be obtained only asymptotically. Additionally, if Vq is not zero then 
\v cannot be constant, since the constraint (2.17) cannot be satisfied in general. Thus, we 
conclude that in general this scenario has Ay ^ 0, and therefore there are not parameter values 
that make it coincide completely with usual (constant mass) massive gravity, as discussed in 
[75]. 

As we observe from Table 1, there exist three critical points and all of them can be 
stable according to the parameter values. Point P\ in the case of standard matter (7 = 1) 
corresponds to a non-accelerating universe, with a dark energy behaving as dust. Although 
it has the advantage that < Qde < 1 ; that is it can alleviate the coincidence problem since 
dark energy and dark matter density parameters can be of the same order, the above features 
disfavor it. Lastly, note that the corresponding graviton mass has become zero. 

Point P2 corresponds to a dark-energy dominated, non-accelerating, universe, with dark- 
energy behaving as dust, and thus it is also disfavored by observations. Additionally, the 
graviton mass remains finite. 

Point P3 is the most interesting solution that can attract the universe at late times. It 
corresponds to a dark-energy dominated universe, which can be accelerating (for < |) or 
non-accelerating according to the parameter values, and where dark energy can lie either in the 
quintessence [87, 88] or in the phantom regime [89] (for Aw Ay < 0). Moreover, the graviton 
mass dynamically becomes zero. These features make this point a very good candidate for 
the description of late-time universe, in agreement with observations. Furthermore, note that 
if the universe starts from the quintessence regime, then the attraction to P3 implies the 
phantom-divide crossing [90] . The realization of the phantom regime and /or of the phantom- 
divide crossing, is a great advantage of extended nonlinear massive gravity, as was analyzed 
in detail in [75, 76]. 

We mention here that naively it looks strange that P3 can be a phantom solution although 
the graviton mass tends to zero and the model should look like quintessence. However, this 
is easily explained since, as we show in Appendix A.l, in the case Aw Ay < where -P3 is 
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Figure 1. Trajectories in the Y-Z plane of the cosmological scenario (3.4), where the varying graviton 
mass square V(ip) is imposed at will, in a flat universe. We use 7 = 1, \y = 2, Aw = —1, «3 = 0:4 = 
0.1. The physical part of the phase space is marked by the shadowed region limited by the red lines. In 
this specific example the universe is led to the phantom stable late-time solution P3. 

phantom, V tends to zero but W and H tend to infinity, which is a Big-Rip-type behavior 
(realized at infinity and not at a finite scale factor) [91-97], that is a typical fate of phantom 
scenarios. In other words, the graviton mass does tend asymptotically to zero, but its previous 
effect has already led the universe to a phantom regime without exit (although not so strongly 
in order to exhibit a Big Rip at a finite scale factor). 

In order to present the above behavior in a more transparent way, we evolve numerically 
the autonomous system (3.4) in the invariant set u = 0, for the the parameters 7 = 1 (dust 
matter), Xy = 2, \\y = — 1, 03 = 04 = 0.1, and in Fig. 1 we depict the corresponding phase- 
space behavior in the Y-Z plane. The physical part of the phase space is marked by the 
shadowed region limited by the red lines. As we observe, in this specific example the universe 
results in the phantom stable late-time solution P3. 

4.1.2 Open universe 

As we show in detail in Appendix A. 2, and as we have depicted in Tables 2, 3 and 4, the 
scenario at hand admits many stable late-time solutions, and this reveals its advantages 
and capabilities, comparing to standard quintessence, as well as to standard (constant-mass) 
nonlinear massive gravity. Note that this scenario admits also curves of solutions apart from 
individual points, which is an additional indication of its generalized features. 

In particular, the first interesting solution, that can attract the universe at late times, 
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is the curve of critical points Q\. It corresponds to an accelerating universe, in which dark 
energy lies always in the phantom regime. Moreover, it has the advantage that < £Ide < 1) 
that is it can alleviate the coincidence problem, and the graviton mass can be zero or not 
according to the parameter values. 

The curves of critical points Q3 and Q^, as well as the individual critical points Q 7 , can 
be stable and thus attract the universe at late times, however since they correspond to zero 
acceleration are not favored by observations (although they have < Qde < 1 an d thus they 
can solve the coincidence problem). 

The curves of critical points Q4 and Q5 can be stable (although their stable manifold 
has smaller dimensionality and thus the stability is weaker) , that is they can be the late-time 
state of the universe, corresponding to an accelerating or non-accelerating universe according 
to the parameter values. Furthermore, note that according to the parameter values they can 
lie in the quintessence or phantom regime, and they possess < Qde < 1- Additionally, 
the graviton mass becomes zero. These features make Q4 and Q5 good candidates for the 
description of the universe. 




0.4 0.6 0.8 

X 



Figure 2. Trajectories in the X-Y plane of the cosmological scenario (3.14) ; where the varying 
graviton mass square V(ip) * s imposed at will, in an open universe. We focus on the invariant set 
&k = U = Z = and we choose 7 = 1, Ay = —2, \\y = 1, 013 = = 0.1. In this specific example the 
stable late-time state of the universe is the phantom solution Qy. 



In particular, as we discussed in paragraph 3.1.2, the curve Q$ contains the quintessence- 
like critical points presented in Table 4, which are obtained in standard quintessence too in 
flat [79] or non-flat geometries [85]. Note however that the stability properties are slightly 
different, since now we have the additional direction of the graviton mass. Amongst these 
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points, Q14 is stable, corresponding to a dark energy-dominated, quintessence universe, which 
can be accelerating or non-accelerating according to the parameter values, and thus it is a 
good candidate for the description of the universe. On the other hand point Q15, which is 
stable in standard quintessence, in the present case it becomes saddle and therefore it cannot 
be the late-time state of the universe. 

Let us present the above results more transparently. In Fig. 2 we show the correspond- 
ing phase-space behavior in the X-Y plane, as it arises from numerical elaboration of the 
autonomous system (3.14). We focus on the invariant set f2/t = U = Z = and we choose 
7 = 1, Ay = — 2, Aw = 1,03 = «4 = 0.1. In this specific example the stable late-time state 
of the universe is the phantom solution Q\. Similarly, in Fig. 3 we depict the corresponding 
phase-space behavior of the autonomous system (3.14), but restricted to the invariant set 
£lk = X = Z = 0, and using 7 = 1, Ay = 2, Xw = 1, «3 = «4 = 0.1 and U c = In this case 
the late-time stable solution of the universe is the quintessence-like point Q14. 




Figure 3. Trajectories in the Y-Z plane of the cosmological scenario (3.14), where the varying 
graviton mass square V(ip) * s imposed at will, in an open universe. We focus on the invariant set 
fife = X = Z = and we choose 7=1, Xy = 2, Xw = 1, 0:3 = 0:4 = 0.1 and U c = In this specific 
example the stable late-time state of the universe is the quintessence-like point Q14. 



Finally, in Fig. 4 we present the phase-space behavior of the autonomous system (3.14), in 
the subset X = Z = 0, which is invariant provided 3 + 3«3 + 04 = 0. In this case the universe 
can be attracted by two stable late-time solutions, namely the expanding, non-accelerating 
Qq (its basin of attraction is the half-subspace > 0), and the contracting Qq (its basin of 
attraction is the half-subspace < 0). Finally, as we discussed in the end of paragraph 3.1.2 
and in Appendix A. 2, we mention that in the scenario at hand the sign of £1^ is invariant. 
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Figure 4. Trajectories of the cosmological scenario (3.14), where the varying graviton mass square 
V(tp) is imposed at will, in an open universe, in the subset X = Z = 0, which is invariant provided 
3 + 3q.3 + Q4 = 0, «3 7^ — 2,Q4 =/= 3. We use the parameters values 7 = 1, Xw = 3. In this specific 
example the stable late-time solutions of the universe are the expanding, non- accelerating Q§ (its 
basin of attraction is the half-subspace Clk > 0), and the contracting Qg (its basin of attraction is 
the half-subspace ttk < 0). Additionally, we can see the saddle points Q15 (non- accelerating with 
< Qde < 1/. Qie (non- accelerating, matter- dominated), Qn (curvature-dominated, contracting) 
and Qi% (non- accelerating, curvature-dominated, expanding), as well as the unstable points Q 12 and 
Q13 (non- accelerating, dark-energy dominated, with stiff wjje)- 

Thus, although our model admits expanding (lower half of Fig. 4) and contracting evolution 
(upper half of Fig. 4), there is no transition from contracting to expanding solutions or vice 
versa, that is a cosmological bounce or turnaround is not possible. 

In summary, as we can see, the scenario of extended nonlinear massive gravity in open 
geometry has a great variety of stable late-time solutions, as was shown in [75, 76] through 
specific examples. 

4.2 Imposing b(t) at will 
4.2.1 Flat universe 

In this case the scenario at hand admits a variety of stable late-time solutions. In particular, 
point R£ corresponds to an expanding dark-energy dominated universe, with dark energy 
lying in the quintessence regime, which can be accelerating or non-accelerating according to 
the usual potential exponent, and the graviton mass is zero. This point exists in standard 
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quintessence too [79], and it is quite important since it possesses wde and q compatible with 
observations. 

Point R§ has the advantage that < Qde < 1) that is it can alleviate the coincidence 
problem, and moreover the graviton mass is zero, however it has the disadvantage that it 
is not accelerating and possesses wde = 0, which are not favored by observations. This 
point exists in standard quintessence too [79], however note that in the present case it is 
non-hyperbolic, and thus its stability is weaker (due to the existence of an extra dimension 
in the phase space, namely the graviton mass). 

Point Rq exists for A^y = and it is always stable. Although at first sight it seems to 
be the \\y —> limit of this is not the case since the complete equations are different. 
It corresponds to an accelerating, dark-energy dominated universe, in which dark energy 
behaves like a cosmological constant, and moreover the graviton mass is zero. 




Figure 5. Trajectories of the cosmological scenario (3. 25) -(3. 28), where the Stuckelberg field function 
function b(t) is imposed at will, in a flat universe, using 7 = 1, Xw = 1jCK3 = ot± = 0.5, B = 1.7. In 
this specific example the stable late-time state of the universe is the expanding, dark-energy dominated, 
quintessence-like point . Additionally, we depict the saddle points R\ (non- accelerating, matter- 
dominated), and R2,Rz (non- accelerating, dark-energy dominated). 



The curves of critical points Rj and Rs can also be the late-time state of the universe (they 
are non- hyperbolic and thus their stability is weaker). They correspond to non-accelerating 
solutions, where the dark energy behaves like dust and where < VLde < 1 5 and additionally 
they possess a non-zero value for the graviton mass. These features disfavor these curves of 
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critical points. Finally, we mention here that although the aforementioned individual points 
were obtained in [86] too, these curves of critical were missed, due to the fact that in the 
analysis one of the phase-space directions was frozen for simplicity. 

In Fig. 5 we depict orbits of the autonomous system (3.25)-(3.28), restricting to the 
invariant set u = 0, and using 7 = l,\w = 1,«3 = «4 = 0.5, B = 1.7. In this specific 
example the stable late-time state of the universe is the expanding, dark-energy dominated, 
quintessence-like point . 

4.2.2 Open universe 

This scenario possesses only one stable solution that can attract the universe at late times, 
namely Sq (although at first sight it seems to be the Xw — > limit of Sf this is not the 
case since the complete equations are different). This point corresponds to a dark-energy 
dominated, accelerating universe, with zero graviton mass, and where the dark energy behaves 
like cosmological constant. This solution is the global attractor of this cosmological system, 
that is the universe will be always led there, for every initial conditions. These features 
make this point a good candidate for the description of the universe. However, we mention 
that it exists only for \\y = 0, that is a form of parameter-tuning is needed. On the other 
hand, for \\y ^ the system does not accept any stable solutions, due to the fact that 




Figure 6. Trajectories of the cosmological scenario (3.34)-(3.38), where the Stiickelberg field function 
function b(t) is imposed at will, in a non-flat universe, restricted to the invariant set u = v = 0, 
using 7=1 and \w =0. In this specific example the stable late-time state of the universe is the 
cosmological- constant-like solution, S§ . Additionally, we depict the saddle points Si (non- accelerating, 
matter- dominated) , and S^S^ (non- accelerating, dark-energy dominated). 
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there are unstable directions related to both curvature and graviton mass. In summary, this 
implies that in general the scenario at hand has disadvantages, unless one tunes the model 
parameters. Finally, note that since the sign of £lk is invariant, although the model admits 
expanding and contracting evolution, a cosmological bounce or a turnaround is not possible. 

In Fig. 6 we present orbits of the autonomous system (3.34)-(3.38), restricting to the 
invariant set u = v = 0, and using 7 = 1, Aw = 0. Note that the evolution is independent 
of the values of 03, 04 and bo, since they do not appear explicitly in the equations governing 
the dynamics in this invariant set. In this specific example the stable late-time state of the 
universe is the cosmological-constant-like solution, Sq. 

5 Conclusions 

In this work we investigated the dynamical behavior of extended (varying-mass) nonlinear 
massive gravity, which is an extension of the usual nonlinear massive gravity [6, 7] where the 
graviton mass is promoted to a scalar-field potential [74]. This scenario has a lot of freedom 
due to the involved free functions, and thus its cosmological implications are significant. 

In order to extract the basic features of the above paradigm, we performed a detailed 
dynamical analysis in the case of an open geometry, adding for completeness the flat case, 
although it proves to have disadvantages that can be cured only at the phenomenological level. 
In both analyses we followed two approaches, namely the theoretically robust one to impose 
the graviton mass square at will and let the equations determine suitably the Stiickelberg field 
function, or the theoretically less-justified one to impose the Stiickelberg field function at will 
and let the equations to determine the graviton mass square. In all cases we extracted the 
late-time solutions and we calculated the corresponding observables, such as the dark-energy 
equation-of-state parameter, the deceleration parameter, and the dark-energy and matter 
density parameter. 

One basic feature of the scenario at hand is that it can lead to an accelerating universe, 
with an effective dark energy lying in the quintessence or in the phantom regime, or experience 
the phantom-divide crossing during the evolution. This is a great advantage since the model at 
hand utilizes only a canonical field. Additionally, and more interestingly, the universe cannot 
only be phantom at one stage of its evolution, but also at its final late-time solutions it can be 
quintessence or phantom like. This is not the case in other modified-gravity scenarios, where 
the universe results to quintessence-like solutions even if it has passed through the phantom 
regime [98]. The above features were discussed in [75, 76] using specific solutions, but in the 
present work they arise from a general dynamical analysis. 

An additional advantage of extended nonlinear massive gravity is that the graviton mass 
goes asymptotically to zero at late times, without fine-tuning, which is in agreement with 
observations. Note that this is not the case in usual massive gravity, where ones needs to 
fine-tune the graviton mass to a very small value by hand. 

Finally, another advantage of the present scenario is that the dark energy density parame- 
ter at the late-time solutions can be between zero and one, which can alleviate the coincidence 
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problem since dark energy and dark matter density parameters can be of the same order. 

In the above analysis we used the exponential ansatz for the usual scalar-field potential, 
and then we used an exponential form for the graviton square mass, in order to be phe- 
nomenologically consistent. One could ask whether the above behaviors are a result of these 
specific ansatzes, or they have a general character. Although this would need an explicit 
investigation from the beginning, the details of our analysis indicate that the results are qual- 
itatively robust for many phenomenologically consistent varying graviton mass choices too. 
However, in the alternative and less-justified approach where the Stiickelberg field is imposed 
at will, our results are quite sensitive to the input ansatz, and therefore a detailed analysis is 
required for every new choice. The fact that the results are very sensitive in the Stiickelberg 
field ansatz, is known to happen in the usual nonlinear massive gravity too [71-73]. 

In summary, the scenario of extended (varying-mass) nonlinear massive gravity, exhibits a 
larger variety and a richer structure of interesting cosmological late-time solutions, comparing 
to usual quintessence, phantom, and quintom cosmology, and also to usual (constant-mass) 
massive gravity. These features are in agreement with observations and thus they make this 
paradigm a good candidate for the description of nature. However, an additional requirement 
for the validity of this scenario is to behave consistently beyond the background level too. 
Since the theory at hand is based on the usual massive gravity formalism in order to become 
Boulware-Deser ghost free, the perturbation analysis could reveal interesting issues too [99, 
100]. Although such a perturbation investigation is therefore necessary, it lies beyond the 
scope of the present work and it is left for a future project. 
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A Stability when V(i(j) is imposed at will 
A.l Flat universe 

For the critical points (u c ,Y c ,Z c ) of the autonomous system (3.4), the coefficients of the 
perturbation equations form a 3 x 3 matrix Q, however since they are quite complicated 
expressions we do not display them explicitly. Despite this complicated form, using the 
specific critical points presented in Table 1, the matrix Q obtains a simple form that allows 
for an easy calculation of its eigenvalues. The corresponding eigenvalues and the stability 
conditions for each critical point are presented in Table 9. 
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Cr. P. 


Eingenvalues 


Stability 


Pi 


-1,3(7-1), 3 (t-tf) 


stable for 7 < min |l, ^ j , X v > § 
saddle point otherwise 


P2 


-1,-3(7-1), 3(l-^) 


stable for 7 > 1, ^ > 1 
saddle point otherwise 


P 3 


-1.-3(7-^), -3(1-^) 


^<min{l,7},A y >f 
saddle point otherwise 



Table 9. The eigenvalues of matrix Q of the perturbation equations of the autonomous system (3.4), 
and the corresponding stability conditions. 



Since in the special case 7 = 1 (dust matter) one eigenvalue of P\ and P2 becomes zero, 
we need to examine this case separately. For 7 = 1 the system (3.4) is restricted to the 
invariant set u = and it admits the general solution 

3 - 2X1 



y( t ) = : v 

eCl( 2A^3)-^I + r_ 2A 



2 
V 



Z(r) = 2 (A.l) 

y/e^ixl+r -2X 2 v e Cl+h ^T 

where c\ and C2 are integration constants. In this case, the system (3.4) admits two classes of 
critical points: the point P3 for which the stability conditions reduce to < 1, X v > | and 
the Z-axis which is a curve of equilibrium points containing the points Pi and P%. The center 
direction of the curve is tangent to the Z-axis, and therefore this curve of critical points is 
normally hyperbolic [101] (a set of non-isolated singular points is called normally hyperbolic 
if the only eigenvalues with zero real parts are those whose corresponding eigenvectors are 
tangent to the set), and since by definition any point on a set of non- isolated singular points 
will have at least one eigenvalue which is zero, all points in the set are non-hyperbolic. The 
stability of a set which is normally hyperbolic can be completely classified by considering the 
signs of the eigenvalues in the remaining directions [101]. In conclusion, in the special case 
7 = 1, the curve of critical points that contains Pi and P2 is stable for > 1. 

Finally, lets us comment on the asymptotic behavior of P3. From the constraint equation 
(2.17) it follows that 

<# _ /i(a re ye~ r ) 
dr Ay/ 2 (a re /e" T )' 

which has the solution 



(A.2) 



Av(V>-<Ao)= [ T { l( " refe Z\ dv, (A.3) 
Jo h(a reJ e 
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where the current scale factor is set to 1 and i/jq denotes the current ^-value. Hence, 

(1 - a re f) [-3a 3 (a re/ - 1) + a 4 (a ref - l) 2 + 3] 



V oc e -VW-^o) 



while 



(e T - a re f) e 2r (3a 3 + 04 + 3) + a^e/ ~~ a re fe T (3as + 2a^) 



(A.4) 



W oc e 



-<W0/>-V>o) 



1 - a re /) [-3a 3 (a re / - 1) + a 4 (a re / - l) 2 + 3] 



e f) e 2r (3a 3 + a 4 + 3) + a 4 a 2 e/ 



e/ e T (3a 3 + 2a 4 ) 



(A.5) 



Therefore, since a re f < 10 9 , if > both V and VF tend to zero as r — > +00. However, 
for < 0, V tends to zero but W tends to infinity as r — > 00, and since Y c / we deduce 
that H — y 00 as t — ^ 00. This is a Big- Rip- type behavior, however it is realized at infinity 
and not at a finite scale factor [91-97]. 

A. 2 Open universe 

Let us discuss the critical points of the autonomous system (3.14) and their stability condi- 
tions. From the last equation of (3.14) it follows that either q = or £1^ = 0, and therefore 
we can simplify the investigation and examine these two cases separately. 

Note that the variable choice (3.13) allows for an easy, partial, classification of expanding 
and contracting solutions. In particular, solutions with Qk = k/{aH) > correspond to 
H > and thus to expansion, while those with 0,^ < correspond to H < and therefore 
to contraction (k = y/\K\ throughout this work). That is why points with > are 
denoted with the subscript "+", while those with < are denoted with the subscript 
However, this is only a partial classification, since it cannot work for solutions with = 0, 
which can be either expanding or contracting. Finally, we mention that although our model 
admits expanding and contracting solutions, since the sign of 0^ is invariant, there is no 
transition from contracting to expanding solutions or vice versa. Nevertheless, there could 
still exist at the edge of the phase space, and in such a case they could be revealed only 
through application of Poincare central projection method [82-84]. This analysis lies beyond 
the scope of the present work and it is left for future investigation. 
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Analysis in the invariant set £lf. = 

In this case, from the first equation of (3.14) it follows that X = 0. Thus, the curvatureless 
equilibrium solutions must satisfy 

Y \VQU 2 [( 7 - 2)Ay - X w ] + 2U(X V X W + 3) 

-V6{X W Y + Ay [7 - 7 Y + (7 - l)Z(4a 3 + a 4 + 6)]}} = 0, (A.6) 
Z {-v^(7 - 3)X V U 2 - 2 (Ay + 3)U 

+VG{X W Y + Ay [7 - 7 y + (7 - l)Z(4a 3 + a 4 + 6)]}} = 0, (A.7) 

U^^-X V . (A.8) 

Note that in this case the evolution equation for U reduces to U' = 0, which implies that in 
the former expressions U behaves as a parameter (a constant). 

Thus, in the case of Sl^ = we have the following curves of critical points: 

. Curve Qf. X cl = 0,Y cl = hvh±h^^±^\ , Zcl = Z c ,U cl = 0,ft fecl = 0, with 
eigenvalues 

{-1,-1,0,0,-37}. 



Curve Q 2 : X c 2 = 0, Y c2 = 0, Z c2 = Z c 

^ + VAt,+6A^[( 7 -3) 7 +(7-3)(7-l 
v / 6(7-3)A v 

' _ 27+A^+VAy+6A^[(7-3)7+(7^3)(7^1)^c(4a 3 +a4+6)+l]+9-3 
U J ' 2(7-3) ' 

\2 



A v + V /A v+ 6A v[(7-3)7+(7-3)(7-l)^ c (4a 3 +a4+6)+l]+9+3 „ , . , 

C/ c2 = — v ^ v / 6(7-3)A v ' fcc2 = ' eigenvalues 



6(7-3) (7- 1) K- Z c (4a 3 +a 4 +6) 



(27-5)A^+^/A4,+6A2,[( 7 -3)7+(7-3)(7-l)Z c (4a 3 +Q:4+6)+l]+9+3 ' 

(A- i /-A w )(A^ + v / A^+6A2 r [(7-3)7+( 7 -3)( 7 -l)Z c (4a 3 +a4+6)+l]+9+3) 

(7-3)Ay 



Curve Q 3 : X c3 = 0, Y c3 = 0, Z c3 = Z c 

^+y / Af / +6A 2 v [( 7 -3)7+(7-3)(7-"i: 
v / 6(7-3)A v 

-27-A^ + v / A^+6A^[(7-3)7+(7-3)(7-l)Zc(4a 3 +a4+6)+l]+9+3 



U c3 = Z^/S^EE^^^0^±^±^Z1 ^ n kc3 = 0, with eigenvalues 



{0,-1, 



2(7-3) 

6(7-3)(7-l)A^Z e (4a 3 +a 4 +6) 



(27-5)A2,-7A^+6A^[(7-3)7+(7-3)(7-l)Z c (4a 3 +Q:4+6)+l]+9+3 ' 

{\v-*w)(->$WK+ 6X V [(7-3)7+(7-3)(7-l)Zc(4a 3 +a4+6)+l]+9-3) 

(7-3)A v J - 

Curve Q4: X C 4 = 0, y c 4 = 0, Z c ^ = 0, U C 4 = U c , Qk C 4 = 0, with eigenvalues 

f n , v / 6(2-37)Av+3v / 6( 7 -2)Ayf/ e 2 +12[/ e 3[-V6 1 \ v +V6(-,-3)X v U^+2(xl+3)U c ] 

3{-v / 67Av+v / 6a c 2 [(7-2)Av-Avy]+2[/ c (A v A X y+3)} \ 
3U C — v/6Ay J" 
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. Curve g 5 : X c5 = 0, Y c5 = 1 - C/ c 2 + (^^zil^Z^) ) Zc5 = 0) u<£ = u& ^ = . 
In order to determine the stability of this curve of critical points we need to resort to 
numerical inspection. 

The examination of the sign of the above eigenvalues is straightforward for the general 
case however in the special case 7 = 1, which is the most interesting in physical terms 

since it corresponds to dust matter, some eigenvalues become zero and thus the corresponding 
curves of critical points become non-hyperbolic. In this case if the remaining eigenvalues have 
different sign then the curve of critical points behaves like saddle, while if they are of the 
same sign then the non-hyperbolic curve of critical points has a stable or unstable manifold of 
smaller dimensionality (in principle one must apply the center manifold theorem [101]). The 
curves of critical points Q1-Q5 for the special case 7 = 1 are summarized in Table 2, while 
their stability conditions are displayed in Table 3. 

Analysis in the invariant set q = 

In the case q = 0, from (3.14) we deduce that the equilibrium solutions must satisfy one of 
the following three possibilities: 

. Y c = 0,Z c ^0,U c = ^, 
• Y c = 0, Z c = 0. 

In the first case, substituting the values of Z c = 0, U c = into the fourth equation of 



(3.14) we conclude that the equilibrium solution satisfies Y c = t^-- Inserting this into the 



3X\y 

3A 



ir 



expression for q we obtain the additional constraint — ^ 7 2 ^ x w^ k — = o, which leads 

to Ofe c = ±. /l — -rf- (corresponding to expanding and contracting universe respectively). 
V w 

Finally, inserting these expressions in the relation for b (3.16) we find that Slkb = X, and thus 
the first equation of (3.14) is satisfied identically, irrespectively the value of X. In summary, 
in this case we obtain two curves of critical points, namely 



Qt '■ — x c , Y c g — ^ 2 , z+q — o, u^q — , r^ c6 — Jl ^ 2 , 



and 



In the second case, the system admits two critical points, namely 



Qt : x+ - o,y+ - o,z+ - -^-^-JL___ £/+ _ ^,n+ c7 - Jl- i 
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and 



V6_ 

3A y (4a 3 + a 4 + 6)'" c7 3Ay 



Ql '■ X c7 - °> Y c7 - °> Z c7 - ~q\2 (a„._ , „. r~cT' ^c7 ~ ^T7T> ^fcc7 




Finally, in the third case, from the fourth equation of (3.14) it follows that U c = 0. Thus, 
substituting Y c = 0, Z c = 0, U c = in the rest of the equations, and assuming that 7 7^ |, we 
obtain Q,k c = ±1. Therefore, for £lk c = +1 the first equation of (3.14) gives 

Q s : X cS = 1, Y c8 = 0, Z c8 = 0, U c8 = 0, O fcc8 = 1 

Q9 '■ X c § = -, Y c g = 0, Z c g = 0, U c g = 0, S7fc c g = 1 

n ■ ^ - V 4a 3- 6a 4 + 2a3 + Q4 _ _ _ 

V10 • A cl0 — ; *cl0 — U, Z c io — U, U c io — U, lifccio — 1, 

«4 



while for ilfc c = —1, we obtain that X = X c n, where X c n is the unique real solution of the 
equation -2a 3 (X 2 + X - 2) + a A {X + 1){X - l) 2 + 6 = 0. 

■ ■ 



Lets us now examine the eigenvalues associated to the critical points or curves 



Qn- The eigenvalues of Q± are (o, 2 - 3 7 , 2 - - V 8A w sx w _ x V 8X ' 2 w sx w _ A The 

eigenvalues of Qf (for 7 = 1) are jo, — 1, 2 — Ai(a3, 04, Ay, Xw), A 2 («3, 04, Ay , Avy)j, 
where Ai^o^, 04, Ay, Aw) are complicated functions of their arguments that can be obtained 
explicitly only by numerical elaboration. The eigenvalues associated to Qs are {0, 2, 2, — 2, 2 — 
37}. The eigenvalues associated to Qg,Qio are {—2,2,2,0^3(7,03,04)} and finally for Qn 
they are {—2, 2, 2, 0, A 4 (7, 03, 04)}, where A3 i 4( 7 , 03, 04) are complicated expressions of their 
arguments. Thus, Qs-Qn are always saddle since at least two eigenvalues have different signs. 

The individual critical points Q7 -Q10 and the curves of critical points Q^ and Qn, for the 
special case 7 = 1, are summarized in Table 2, while their stability conditions are displayed 
in Table 3. 

Quintessence-like solutions 

We close this Appendix by mentioning that the curve of critical points Q5 analyzed above 
includes many interesting cosmological solutions, and in particular the points of standard 
quintessence [79, 85]. Focusing for simplicity on the case 7 = 1, these points were presented 
in Table 4. However, the stability conditions are different than the usual conditions in [79, 85] 
due to the presence of extra phase-space directions, namely those of curvature and graviton 
mass. 

In particular, for the critical points Q12 to Qis of Table 4, the coefficients of the pertur- 
bation equations form a 5 x 5 matrix Q, that allows for an easy calculation of its eigenvalues. 
The corresponding eigenvalues and the stability conditions for each critical point are dis- 
played in Table 10. Finally, some of these points possess one zero eigenvalue and are thus 
non-hyperbolic. In the case of normally-hyperbolic curves of critical points (that is the only 
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Cr. P. 


Eingenvalues 


Stability 


Ql2 


2, -1, 0, y/6X v + 6, \/~6X w + 6 


saddle point 


Ql3 


2, —1, 0, 6 — \/6Ay, 6 — \/6Xw 


saddle point 


Qu 


1'°' 2\ v -\ w ' 2 (. A W 2 J,A^(A W Ay) 


stable node for 
— V 2 < Avk < 0, Ay < Xw or 
U < Aw < V Z, Ay > Aiy 
saddle point otherwise 


Ql5 


i 1 n 3(Ay — Avv) 9(A\/ — A;y) 
X > 2' U ' X w ' A vv (2AvA w -3) 


111 j_ 

saddle point 


Ql6 


3,3,-l,|,0 


saddle point 


Ql7 


2,2,-1,0,2 


saddle point 


Ql8 


2,2,-1,-1,-4 


saddle point 



Table 10. The eigenvalues of matrix Q of the perturbation equations of the autonomous system (3.14), 
and the corresponding stability conditions, for the quintessence-like solutions presented in Table 4. 

eigenvalues with zero real parts are those whose corresponding eigenvectors are tangent to the 
set) the stability is extracted considering the signs of the rest eigenvalues [101]. For isolated 
non-hyperbolic critical points we can determine the dimensionality of their stable manifold 
using the linearization technique [101]. 

B Stability when b(t) is imposed at will 
B.l Flat universe 

For the critical points (x c , y c , u c , v c ) of the autonomous system system (3.25)-(3.28), the coef- 
ficients of the perturbation equations form a 4 x 4 matrix Q, which using the specific critical 
points presented in Table 5 acquires a simple form that allows for an easy calculation of its 
eigenvalues. The corresponding eigenvalues and the stability conditions for each critical point 
are presented in Table 11. We mention that point R± is non- hyperbolic, but since it has 
eigenvalues with different sign, and using the center manifold analysis [101], we can straight- 
forwardly show that it behaves as saddle point. Moreover, the non-hyperbolic curve of critical 
points Rj has a central direction normal to the set and therefore it behaves as stable. Finally, 
note that although R§ at first sight seems to be the Xw — > limit of Rf this is not the case 
since the complete equations are different. 
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Cr. P. 


Eigenvalues 


Stability 


R\ 


_3 3 1 n 

2 ' 2 ' ' 


non-hyperbolic (behaves as saddle point) 


R2 


3, 3, —1, 3 — y^Xw 


saddle point 


R3 


3, 3, —1, 3 + </ |Avy 


saddle point 


Ra 


^ 2 2 
— 1, — 3 H — q , — 3 + Ati/, — 3 + A11/ 


stable node for Aw < 3 
w 

saddle point for 3 < X w < 6 


Rf 


0,-1, a~ (Aw), a + (Xw) 


non- hyperbolic 
3D stable manifold for 
3 < X w < f or 

\2 - 24 
A W > T 




-3,-3,-3,-1 


stable node 


R 7 


-1, 0, /3~(\ w , He), /3 + (Aiy, y c ) 


normally hyperbolic (behaves as stable) 


Rs 


—3, —1, 0, | — ^J\\wx c 


stable for x c Xw > y^|; 
saddle point otherwise. 



Table 11. The eigenvalues of matrix Q of the perturbation equations of the autonomous sys- 
tem (3.25)-(3.28), calculated at the critical points presented in Table 5, and their stability condi- 

and ^{XwiVc) = 



tions. We have introduced the notations a ± (Avi/) 
-| [3- A 2 w yl ± V^-yc 4 -18(A^-2)^ + 9 



-1 ± 



B.2 Open universe 

For the critical points (x c ,yc,Uc,v c ,^kc) of the autonomous system system (3.34)-(3.38), the 
coefficients of the perturbation equations form a 5 x 5 matrix Q, which using the specific 
critical points presented in Table 7 acquires a simple form that allows for an easy calculation 
of its eigenvalues. The corresponding eigenvalues for each critical point are presented in Table 
12. Finally, note that although at first sight seems to be the Xw limit of S^ 1 this is 
not the case since the complete equations are different. 

In order to examine the corresponding stability conditions we have to examine the sign 
of these eigenvalues. An interesting observation from (3.37) is that the sign of v (which 
according to (3.33) is the auxiliary variable proportional to the graviton mass square) is 
invariant. Therefore, v remains zero if initially it is zero, and in this case we can examine the 
system in the invariant set v = 0. In this case the possible late-time solutions are either 5^ 
provided X w < 2 or either for X w > 2. In the particular case of 2 < A^ < |, the points 

are spiral attractors in a 2D sub-manifold (two negative real eigenvalues and two complex 
conjugated eigenvalues with negative real part). Finally, points S^ 1 are non-hyperbolic, with 
a 4D stable manifold. 
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Cr. P. 


Eigenvalues 


Stability 


Si 


-11 -1 3 1 

2 ' 2 ' ' ' 2 


saddlp noint 


S2 


d 10 /3\ 
0, Z, O, —1, — \ 2 W 


saddle point 


s 3 


6, 2, 3, —1, 3+ -v/|Aiy 


saddle point 


°4 


\2 1 1 A w 1 \2 1 1 A w 
A W> ' ~r 2 ' ' A W> ' 2 


saddle point 


q± 


-l,a-(X w ),a + (X w ),3, \ 


saddle point 


st 


0,-3,-3,-3,-1 


non-hyperbolic (4D stable manifold) 


sf 


-2,2,-1,-1,1 


saddle point 


S 8 




saddle point 



Table 12. The eigenvalues of matrix Q of the perturbation equations of the autonomous system 
(3.34)-(3.38), calculated at the critical points presented in Table 7, and their stability conditions. We 

/ /24A 2 - — 7A 4 \ 

have introduced the notations q ± (Ah?) = | I — 1 ± — p — J . 



However, in the case where v 7^ only points S^ behave as stable, since all the other be- 
come saddle points. In particular, introducing the local coordinates {x — x c ,y — y c , u, v, ilk} = 
e jx, y, u, v, f2fc j + 0(e) 2 where e is a constant satisfying e < 1, we deduce that 

v' = Zv (x 2 - y 2 c + 1) + h.o.t 

n k ' = (3x 2 - 3y 2 + 1) - [/? 2 (a 3 + a 4 ) - 2/3(2a 3 + a 4 + 1) + 3a 3 + a 4 + 3] + h.o.t, 

(B.l) 

where x c and y c are the coordinates of the critical points Si to 5*5 and h.o.t denoting "higher 
order terms". These equations admit the general solutions 



v = cqe 



3r(x 2 -y 2 c + l) 



fin 



e 2 



ir(3x2- 3 y2+l) _ 3r(x 2 -y 2 +l) 



(a 3 /3 2 - 4a 3 /3 + 3a 3 + a 4 /3 2 - 2a 4 /3 + a 4 - 2/3 + 3) 



3x 2 — 3y 2 + 5 



+ c 2 e2 



±r(3x 2 c -3y 2 +l) 



(B.2) 



which implies that the system is unstable in v and Qk directions. 

In the special case of point Sq, using a similar approach we extract that the perturbations 
v and 0^ satisfy the equations 



$Jj2 , 1 _ 

fit. ^ 



(B.3) 
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where •& = f3 [/3 2 (a3 + a^) — 2(3(2as + 04 + 1) + 3a^ + + 3] . The system (B.3) admits two 
general solutions 

4e 2c 2 2i?e 2c2 - T 

V ~ ( e r_ e c 2ci1? )2' e r_ eC2c ^ ^ 

and 

4e 2c2 ^ = 2e^ 



(e c 2ci?? + e T ) 2 ' (e c 2cii? + e T ) (2e c 2cii9 + e r 

where ci and C2 are integration constants. In both cases the u-perturbations and f2fc- 
perturbations decay to zero in the limit r — > +00, and thus points are stable. 

The stability conditions for the critical points Si-S^ are summarized in Table 12. 
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